#include "refloat.h"
#include "cxfloat.h"
#include <sstream>
#include <iostream>
#include <stdlib.h>
#include <algorithm>

using namespace std;


thread_local refloat *refloat_cache1=NULL;
thread_local bool refloat_all_allocations_temporary=false;
thread_local std::vector<__mpfr_struct*>* temporaries=NULL;


refloat::~refloat()
{
  if((is_temp==false) || (!temporaries)){free(value->_mpfr_d);free( value );}
  else
  {
    temporaries->push_back(value);
  }
}


void refloat::set_temporary(mpfr_prec_t prec)
{
  if(temporaries->size() == 0)
  {
    allocate_mpfr_struct(value);
    init_mpfr_value(value, prec);
  }
  else
  {
    value = temporaries->back();
    if(prec != mpfr_get_prec(value))
    {
      free(value->_mpfr_d);free( value );
      allocate_mpfr_struct(value);
      init_mpfr_value(value, prec);
    }
    temporaries->pop_back();
  }
}


void init_refloat_cache(mpfr_prec_t p)
{
  bool oldtemp = refloat_all_allocations_temporary;
  refloat_all_allocations_temporary = false;
  refloat_cache1 = new refloat(0,p);
  refloat_all_allocations_temporary = oldtemp;
  
  temporaries = new vector<__mpfr_struct*>();
}


void set_refloat_precision(mpfr_prec_t p)
{
  refloat_cache1->set(0,p);
}


void clear_refloat_cache()
{
  if(refloat_cache1 != NULL)
  {
    delete refloat_cache1;
  }
  
  if(temporaries != NULL)
  {
    for(unsigned long int i=0;i<temporaries->size();++i)
    {
      free((*temporaries)[i]->_mpfr_d);
      free( (*temporaries)[i] );
    }
    temporaries->clear();
    delete temporaries;
  }
  
  
  refloat_cache1=NULL;
  
  temporaries=NULL;
}


void init_refloat_cache(void *arg)
{
  mpfr_prec_t p = *((mpfr_prec_t*)arg);
  init_refloat_cache(p);
}


void set_refloat_precision(void *arg)
{
  mpfr_prec_t p = *((mpfr_prec_t*)arg);
  refloat_cache1->set(0,p);
}


void clear_refloat_cache(void *arg)
{
  clear_refloat_cache();
}


void refloat_finishUp(void *arg)
{
  clear_refloat_cache();
  mpfr_free_cache();
}


void refloat_startUp(void* arg)
{
  init_refloat_cache(arg);
}


void refloat::print(string &re, int base, mpfr_prec_t prec, bool bold) const
{
  stringstream ss;
  ss.clear();ss.str("");
  string s;
  
  char *tstr = new char[(int)(mpfr_get_prec(value))+2];
  mp_exp_t *expt = new mp_exp_t;
  mpfr_get_str(tstr, expt, base, prec, value, GMP_RNDN);
  
  if(bold == true)
  {
    if(tstr[0]=='-')
    {
      s=tstr;s.erase(0,1);
      ss<<"-"<<s[0]<<".";
      s.erase(0,1);
      if(base!=10){ss<<s<<"\033[1m@"<<((*expt)-1)<<"\033[0m";}
      else{ss<<s<<"\033[1mE"<<((*expt)-1)<<"\033[0m";}
    }
    else
    {
      s=tstr;
      ss<<s[0]<<".";
      s.erase(0,1);
      if(base!=10){ss<<s<<"\033[1m@"<<((*expt)-1)<<"\033[0m";}
      else{ss<<s<<"\033[1mE"<<((*expt)-1)<<"\033[0m";}
    }
  }
  else
  {
    if(tstr[0]=='-')
    {
      s=tstr;s.erase(0,1);
      ss<<"-"<<s[0]<<".";
      s.erase(0,1);
      if(base!=10){ss<<s<<"@"<<((*expt)-1);}
      else{ss<<s<<"E"<<((*expt)-1);}
    }
    else
    {
      s=tstr;
      ss<<s[0]<<".";
      s.erase(0,1);
      if(base!=10){ss<<s<<"@"<<((*expt)-1);}
      else{ss<<s<<"E"<<((*expt)-1);}
    }
  }
  
  re=ss.str();
  
  delete[] tstr;
  delete expt;
}


void Si_Ci_Acton(const refloat& z, refloat& Si, refloat& Ci, unsigned long int n)
{
  mpfr_prec_t prec = mpfr_get_prec(Si.value);
  refloat counter(0, prec);
  counter.set(n);
  
  refloat A(1,prec);
  refloat B(0,prec);
  refloat C(0,prec);
  
  refloat temp1(0,prec);
  refloat temp2(0,prec);
  
  while(mpfr_cmp_si(counter.value, 0) > 0)
  {
    temp1.set(counter);
    temp1.mul_power_2(1);
    temp1.reciprocal();
    temp1.mul(z);
    temp1.mul(C);
    A.add(temp1);
    
    B.mul(counter);
    B.mul_power_2(1);
    mul(A, z, temp1);
    B.add(temp1);
    temp1.set(counter);
    temp1.mul_power_2(1);
    temp1.sub(1);
    B.div(temp1);
    
    C.add(B);
    
    temp1.set(C);
    temp1.reciprocal();
    A.mul(temp1);
    B.mul(temp1);
    C.set(1);
    
    counter.sub(1);
  }
  
  mul(z, C, temp1);
  temp1.reciprocal();
  A.mul(temp1);
  B.mul(temp1);
  
  
  
  sin_cos(z, temp1, temp2);
  
  mul(B, temp1, Ci);
  mul(A, temp2, C);
  Ci.sub(C);
  
  mul(B, temp2, Si);
  mul(A, temp1, C);
  Si.add(C);
  Si.neg();
  mpfr_const_pi(temp1.value, GMP_RNDN);
  temp1.mul_power_2(-1);
  Si.add(temp1);
}


void Si_Ci_Acton(const refloat& z, refloat& Si, refloat& Ci)
{
  mpfr_prec_t prec = mpfr_get_prec(Si.value);
  
  refloat tol(1,prec);
  tol.mul_power_2(-((long int)prec));
  
  refloat prev_si(0, prec);
  refloat prev_ci(0, prec);
  refloat diff(0, prec);
  
  unsigned long int n = 16;
  
  while(true)
  {
    Si_Ci_Acton(z, Si, Ci, n);
    
    bool si_conv = false;
    bool ci_conv = false;
    
    sub(Si, prev_si, diff);
    diff.div(Si);
    if(compare_abs(diff, tol) < 0){si_conv = true;}
    sub(Ci, prev_ci, diff);
    diff.div(Ci);
    if(compare_abs(diff, tol) < 0){ci_conv = true;}
    if((si_conv == true) && (ci_conv == true)){break;}
    
    n *= 2;
    prev_si.set(Si);
    prev_ci.set(Ci);
  }
}


void Si_Ci_Continued_Fraction(const refloat& z, refloat& Si, refloat& Ci)
{
  //continued fraction :
  //[z; 1,z,1/2,z,1/3,z...]
  
  mpfr_prec_t target_prec = mpfr_get_prec(Si.value);
  mpfr_prec_t oldprec = mpfr_get_prec(cxfloat_cache1->real_part);
  mpfr_prec_t prec = (3*target_prec)/2;
  bool prec_changed=false;
  if((3*target_prec)/2 > oldprec){set_cxfloat_precision(prec);prec_changed=true;}
  
  cxfloat tol(1,0,prec);
  tol.mul_power_2(-((long int)target_prec));
  
  refloat retemp1(0,prec);
  cxfloat iz(0,0,prec);
  iz.set(retemp1, z);
  
  cxfloat p1(0,0,prec);
  cxfloat p2(1,0,prec);
  cxfloat p3(0,0,prec);
  cxfloat q1(1,0,prec);
  cxfloat q2(0,0,prec);
  cxfloat q3(0,0,prec);
  
  cxfloat* p_current = &p3;
  cxfloat* p_prev = &p2;
  cxfloat* p_prev2 = &p1;
  cxfloat* q_current = &q3;
  cxfloat* q_prev = &q2;
  cxfloat* q_prev2 = &q1;
  
  
  unsigned long int count = 1;
  cxfloat temp1(0,0,prec);
  cxfloat running(1,0,prec);
  cxfloat diff(0,0,prec);
  
  unsigned long int niter = 4;
  
  
  //first calculate E_1(i*z)
  while(true)
  {
    for(unsigned long int i=0;i<niter;i++)
    {
      if((count%2) == 1){temp1.set(iz);}
      else
      {
        temp1.set(count>>1, 0);
        temp1.reciprocal();
      }
      mul(*p_prev, temp1, *p_current);
      p_current->add(*p_prev2);
      mul(*q_prev, temp1, *q_current);
      q_current->add(*q_prev2);
      
      swap(p_current, p_prev2);
      swap(p_prev2, p_prev);
      swap(q_current, q_prev2);
      swap(q_prev2, q_prev);
      count++;
    }
    div(*p_prev, *q_prev, temp1);
    sub(temp1, running, diff);
    diff.div(running);
    magnitude(diff, running);
    if(compare_real(running, tol) < 0){running.set(temp1);break;}
    running.set(temp1);
    niter *= 2;
  }
  cxfloat E1iz(0,0,prec);
  temp1.set(iz);
  temp1.neg();
  exp(temp1, E1iz);
  E1iz.div(running);
  
  
  //now calculate E_1(-iz)
  iz.neg();
  niter=4;
  count=1;
  running.set(1,0);
  diff.set(0,0);
  temp1.set(0,0);
  p1.set(0,0);
  p2.set(1,0);
  p3.set(0,0);
  q1.set(1,0);
  q2.set(0,0);
  q3.set(0,0);
  
  p_current = &p3;
  p_prev = &p2;
  p_prev2 = &p1;
  q_current = &q3;
  q_prev = &q2;
  q_prev2 = &q1;
  while(true)
  {
    for(unsigned long int i=0;i<niter;i++)
    {
      if((count%2) == 1){temp1.set(iz);}
      else
      {
        temp1.set(count>>1, 0);
        temp1.reciprocal();
      }
      mul(*p_prev, temp1, *p_current);
      p_current->add(*p_prev2);
      mul(*q_prev, temp1, *q_current);
      q_current->add(*q_prev2);
      
      swap(p_current, p_prev2);
      swap(p_prev2, p_prev);
      swap(q_current, q_prev2);
      swap(q_prev2, q_prev);
      count++;
    }
    div(*p_prev, *q_prev, temp1);
    sub(temp1, running, diff);
    diff.div(running);
    magnitude(diff, running);
    if(compare_real(running, tol) < 0){running.set(temp1);break;}
    running.set(temp1);
    niter *= 2;
  }
  cxfloat E1niz(0,0,prec);
  temp1.set(iz);
  temp1.neg();
  exp(temp1, E1niz);
  E1niz.div(running);
  
  sub(E1niz, E1iz, temp1);
  temp1.mul_by_i();
  temp1.mul_power_2(-1);
  temp1.get_real(Si);
  mpfr_const_pi(retemp1.value, GMP_RNDN);
  retemp1.mul_power_2(-1);
  Si.add(retemp1);
  
  add(E1iz, E1niz, temp1);
  temp1.get_real(Ci);
  Ci.mul_power_2(-1);
  Ci.neg();
  
  if(prec_changed==true){set_cxfloat_precision(oldprec);}
}


void Si_Ci_power_series(const refloat& z, refloat& Si, refloat& Ci)
{
  mpfr_prec_t target_prec = mpfr_get_prec(Si.value);
  mpfr_prec_t oldprec = mpfr_get_prec(cxfloat_cache1->real_part);
  mpfr_prec_t prec = (3*target_prec)/2;
  bool prec_changed=false;
  if((3*target_prec)/2 > oldprec){set_cxfloat_precision(prec);prec_changed=true;}
  
  refloat tol(1, prec);
  tol.mul_power_2(-((long int)target_prec) - 2);
  
  Si.set(0, prec);
  Ci.set(0, prec);
  
  refloat temp1(0, prec);
  refloat temp2(0, prec);
  refloat temp3(0, prec);
  refloat temp4(0, prec);
  refloat z2(0, prec);
  z2.set(z);z2.sqr();
  refloat diff(0, prec);
  long int sign_val = 1;
  
  unsigned long int n = 0;
  
  //sine integral
  //temp1 is what we add to Si
  temp1.set(z);
  //temp4 keeps track of the power of z
  temp4.set(z);
  //temp2 keeps track of the factorial
  temp2.set(1);
  while(true)
  {
    Si.add(temp1);
    
    if( compare_abs(temp1, tol) < 0 ){break;}
    
    temp3.set(2*n + 2);
    temp2.mul(temp3);
    temp3.set(2*n + 3);
    temp2.mul(temp3);
    temp4.mul(z2);
    temp3.mul(temp2);
    sign_val = -sign_val;
    
    temp1.set(sign_val);
    temp1.mul(temp4);
    temp1.div(temp3);
    
    n++;
  }
  
  temp4.set(z2);
  temp2.set(2);
  temp3.set(4);
  sign_val = -1;
  temp1.set(sign_val);
  temp1.mul(temp4);
  temp1.div(temp3);
  n = 1;
  while(true)
  {
    Ci.add(temp1);
    
    if( compare_abs(temp1, tol) < 0 ){break;}
    
    temp3.set(2*n + 1);
    temp2.mul(temp3);
    temp3.set(2*n + 2);
    temp2.mul(temp3);
    temp4.mul(z2);
    temp3.mul(temp2);
    sign_val = -sign_val;
    
    temp1.set(sign_val);
    temp1.mul(temp4);
    temp1.div(temp3);
    
    n++;
  }
  
  mpfr_const_euler(temp1.value, GMP_RNDN);
  Ci.add(temp1);
  temp1.set(z);log(temp1, temp2);
  Ci.add(temp2);
  
  
  Si.set_prec(target_prec);
  Ci.set_prec(target_prec);
  if(prec_changed==true){set_cxfloat_precision(oldprec);}
}


void Si_Ci(const refloat& z, refloat& Si, refloat& Ci)
{
  mpfr_prec_t prec = mpfr_get_prec(Si.value);
  refloat temp1(0, prec);
  refloat temp2(0, prec);
  
  bool sign_changed = false;
  refloat zz(0,prec);zz.set(z);
  if(compare(zz, temp1) < 0)
  {
    sign_changed = true;
    zz.neg();
  }
  
  temp1.set((unsigned long int)(prec));
  temp2.set("0.5", 10);
  temp1.mul(temp2);
  temp2.set(20);
  temp1.add(temp2);
  if(compare_abs(z, temp1) < 0)
  {
    if(sign_changed == false){return Si_Ci_power_series(zz, Si, Ci);}
    else
    {
      Si_Ci_power_series(zz, Si, Ci);
      Si.neg();
      return;
    }
  }
  
//   temp1.set((unsigned long int)(prec));
//   temp2.set("2.5", 10);
//   temp1.mul(temp2);
//   temp2.set(10);
//   temp1.add(temp2);
//   if(compare_abs(z, temp1) < 0)
//   {
//     if(sign_changed == false){return Si_Ci_Continued_Fraction(zz, Si, Ci);}
//     else
//     {
//       Si_Ci_Continued_Fraction(zz, Si, Ci);
//       Si.neg();
//       return;
//     }
//   }
  
  if(sign_changed == false){return Si_Ci_Acton(zz, Si, Ci);;}
  else
  {
    Si_Ci_Acton(zz, Si, Ci);
    Si.neg();
    return;
  }
}









